function elements = Initialize(elements,model)
ne = model.countElements();
dim = 2;
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    nn = f.countModes();
    ng = e.getNumGaussPoint();
    C = elements(i).constitutiveMatrix;
    for j=1:ng
        xi = e.getGaussPoint(j-1);
        w = e.getGaussWeight(j-1);
        jac = f.jacXAt(xi);
        detJac = det(jac);
        G = computeEnhancedStrainOperator(e,xi);
        B = zeros(3,dim*nn);
        B = e.computeB(f,B,xi);
        if j==1
            H = w*G'*C*G*abs(detJac);
            Gamma = w*G'*C*B*abs(detJac);
        else
            H = H + w*G'*C*G*abs(detJac);
            Gamma = Gamma + w*G'*C*B*abs(detJac);
        end
    end
    elements(i).H = H;
    elements(i).Gamma = Gamma;
end
end
